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For a thermal field theory formulated in the grand canonical ensemble, the distribution of the 
total momentum is an observable characterizing the thermal state. We show that its cumulants 
are related to thermodynamic potentials. In a relativistic system for instance, the thermal variance 
of the total momentum is a direct measure of the enthalpy. We relate the generating function of 
the cumulants to the ratio of (a) a partition function expressed as a Matsubara path integral with 
shifted boundary conditions in the compact direction, and (b) the ordinary partition function. In 
this form the generating function is well suited for Monte-Carlo evaluation, and the cumulants can 
be extracted straightforwardly. We test the method in the SU(3) Yang-Mills theory and obtain the 
entropy density at three difi'erent temperatures. 
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I. INTRODUCTION 

Thermal field theory is a theoretical tool of central im- 
portance in condensed matter physics, plasma physics, 
nuclear physics and cosmology [3, Obtaining first- 
principles predictions from a thermal field theory is of- 
ten challenging, since it describes an infinite number of 
degrees of freedom subject to both quantum and ther- 
mal fluctuations. New theoretical concepts and more ef- 
ficient computational techniques are still needed in many 
contexts, particularly when weak-coupling methods are 
inapplicable. 

In this work we exploit the global symmetries of a ther- 
mal theory to define the contributions to the partition 
function of states with given quantum numbers, and show 
that they are well suited for ah initio Monte-Carlo com- 
putations. Our main results, which are based on spatial 
translation invariance, concern the relative contribution 
to the partition function of states with total momentum 
p. These contributions form the probability distribution 
of p, whose cumulants can be related in a simple manner 
to the equation of state by using Lorentzian or Galilean 
invariance. For instance, in Quantum Chromodynamics 
(QCD) at zero baryon chemical potential ^, the variance 
of the total momentum measures the entropy of the sys- 
tem. As we will see shortly, the generating function K of 
the cumulants of the momentum distribution can be ex- 
pressed as a ratio of two partition functions. The latter 
are represented by Euclidean path integrals with different 
boundary conditions for the fields in the time direction, 
and their ratio can therefore be computed by standard 
Monte Carlo techniques. 

While the following theoretical discussion extends to a 
wide class of theories, we illustrate our ideas numerically 
in the SU (3) Yang-Mills theory, where we determine the 
entropy density of the system at three different temper- 
atures. Crucially, the method can be applied to QCD, 
since the hermiticity properties of finite-difference opera- 



tors are not affected by our 'shifted' boundary conditions, 
and hence the actions to be used in the simulations are 
always real at /z = 0. As an example of an application 
to a non-relativistic system, it can also be employed in 
the study of neutron matter in the Euclidean approach 
of %\. 

There are already several established methods to com- 
ute the thermodynamic properties of gauge theories [3- 
. They require either a vacuum subtraction or a renor- 
malization constant to be determined, facts which make 
it difficult to apply them at arbitrarily high and low tem- 
peratures. Our method avoids these problems, and has 
the further advantage that a Symanzik improvement of 
the action (i.e. a suppression of the discretization er- 
rors @, II]) automatically leads to a corresponding im- 
provement in the thermodynamic quantities. Computa- 
tionally, the method is rather expensive, since it consists 
in calculating a ratio of partition functions. However 
our experience shows that, even with a simple algorithm, 
physics results can be obtained with commonly available 
computing resources. 



II. MOMENTUM DISTRIBUTION 

Recent progress in lattice field theory has made it pos- 
sible to define and compute by Monte Carlo simulations 
the relative contribution to the partition function due to 
states carrying a given set of quantum numbers associ- 
ated with exact symmetries of a theory ilO^, JJ|. Here 
we apply these techniques to a finite temperature and 
density system in the grand canonical ensemble (or the 
canonical ensemble if there is no conserved particle num- 
ber). The relative contribution to the partition function 
of the states with momentum p is given by 
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where the trace is over all the states of the Hilbert space, 
P^P') is the projector onto those states with total mo- 
mentum p, (3 = l/T is the inverse temperature, H the 
Hamiltonian, N the particle number and L the linear 
size of the system. The generating function K associated 
with the momentum distribution is defined as 
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The connected cumulants are obtained from it as follows. 
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where they have been normalized so as to have a finite 
limit when L — > oo and the /i) dependence has been 
suppressed. Finally, we can write 
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where Z{|3,^i,z) = j'^ {(,-f){H-^.N)+ipzy -g partition 
function in which states of momentum p are weighted by 
a phase e''^'^. The shifted partition function Z{/3,iJ.,z) 
can be expressed as a path integral in Euclidean time by 
adopting the 'shifted' boundary conditions 



a;) = ±0(0, a; + 2;) 



(5) 



with the +(— ) sign for bosonic (fermionic) fields respec- 
tively. From Eqs. ^ and ([U, the generating function 
can be written as the ratio of two partition functions. 
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i.e. two path integrals with the same action but differ- 
ent boundary conditions. In a renormalizeable theory it 
therefore has a finite and universal continuum limit. 

We remark that in the large volume regime, the mo- 
mentum distribution can be expressed at each value of 
p as a saddle point expansion. Its leading term in 
for a system with a finite correlation length, is a Gaus- 
sian with a width equal to the second cumulant. Near a 
second-order phase transition on the other hand, the rel- 
ative size of the fourth and second cumulants can serve 
to characterize a universality class. 



III. CONNECTION TO THERMODYNAMICS 



m of the particles) coincides with the particle number 
flux fil. The relation 



/d3a;e^'=3-3 (^3(^) ^3(0)) 

_/ Jd^xe'''^^^Tooix)T33iO))c (rel.) 

m J (Pxe^''-^''^ (n(a;)T33(0))c (non-rel.) 
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then follows from the Ward identities (Wis) associated 
with (a) the conservation of momentum and (b) the con- 
servation of energy (particle number) respectively in the 
relativistic (non-relativistic) case. Eq. (O is valid for any 
xo ^ and non- vanishing momentum in a periodic box 
of size L, and we now take the limit L — 00, with ky,L 
fixed. In a fluid at rest the pressure p{T) is given by the 
thermal average of the stress tensor {Tij) — Sijp{T), and 
the relations ^ (T, ji) = n, (T, /i) — s, Ts + iin — e+p 
hold in the thermodynamic limit. We thus obtain 
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Alternatively, Eq. ([5]) can be derived directly at — 
from the same family of Wis. In this way an expression 
for the finite- volume correction is obtained [Tsj . 

Relativistic theories at /i — constitute an important 
special case that we investigate in more detail, since it 
is relevant both to heavy-ion collisions and to the 
physics of the early universe [l^ [l^. First, the relation 
e + p = Ts implies that the thermal variance of the mo- 
mentum is a direct measure of the entropy of the system. 
Secondly, by establishing a recursion relation among the 
cumulants based on the Ward identities, it is possible to 
show that the fourth order cumulants are related to the 
specific heat of the system [l^, 
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Combining second and fourth cumulants, one can thus 
obtain the speed of sound — One may prove [l^ 
that in conformal field theories, tlie generating function 
is completely determined by its first non-zero cumulant. 
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Thirdly, it can be shown [l3j that finite volume effects in 
^2,0,0 are exponentially small in m(T) ■ L, where m(T) is 
the lightest screening mass in the theory, and the prefac- 
tor is explicitly known fl7'|. In summary, the cumulants 
can be used to calculate thermodynamic properties. 



In thermal field theories, there is a connection between 
the cumulants of the momentum distribution and ther- 
modynamic functions. In relativistic theories, the reason 
is that the momentum density tt^ = T^i can be chosen 
to coincide with the energy fiux. In non-relativistic theo- 
ries, the particles are the only carriers of momentum, and 
therefore the momentum density (divided by the mass 



NUMERICAL IMPLEMENTATION AND 
RESULTS 
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There are many potential applications of formulae ()8l9p . 
As an illustration we consider the SU{3) lattice Yang- 
Mills theory defined in a cubic box of volume with 
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Lat /3/a L/a ro/a K{p.,z.,a) 



Ai 


5.9 


4 


12 


2 


4.48(5) 


17.20(11) 


5.10(3) 


Ala 


5.9 


4 


16 


2 


4.48(5) 


40.71(15) 


5.089(19) 


A2 


6.024 


5 


16 


2 


5.58(6) 


13.05(10) 


4.98(4) 


A3 


6.137 


6 


18 


3 


6.69(7) 


7.32(8) 


4.88(6) 


A4 


6.337 


8 


24 


4 


8.96(9) 


4.32(16) 


5.12(19) 


As 


6.507 


10 


30 


5 


11.29(11) 


2.62(17) 


4.9(3) 


Bi 


6.572 


4 


12 


2 


12.28(12) 


22.22(11) 


6.58(3) 


Bla 


6.572 


4 


16 


2 


12.28(12) 


53.47(16) 


6.684(20) 


B2 


6.747 


5 


16 


2 


15.34(15) 


17.11(15) 


6.53(6) 


B3 


6.883 


6 


18 


3 


18.14(18) 


9.61(9) 


6.40(6) 


B^ 


7.135 


g 


24 


4 


24.5(3) 


5.42(17) 


6 42('20'l 


Bs 


7.325 


10 


30 


5 


30.7(4) 


3.32(18) 


6.1(3) 


Ci 


7.234 


4 


16 


4 


27.6(3) 


57.44(25) 


7.18(3) 


C2 


7.426 


5 


20 


5 


34.5(4) 


36.5(4) 


7.13(8) 


C3 


7.584 


6 


24 


4 


41.4(5) 


24.7(4) 


6.94(12) 



TABLE I. Lattice parameters and numerical results with 
z = (2, 0, 0). The quantity in the last column is the estimator 
for the entropy density, see Eq. (I15p . 



periodic boundary conditions, at temperature /? = 1/T 
and lattice spacing a. Our goal is to show that the en- 
tropy density can be computed in the thermodynamic 
and continuum limit. 

Even though the lattice breaks continuous translation 
invariance, discrete translations remain as exact symme- 
tries of the lattice action. They suffice to define K at 
finite lattice spacing precisely as in section dlTl. The 
function K{f3, z, a) defined on the lattice differs by O(a^) 
effects from its continuum counterpart. 

The canonical partition function is given as usual by 

Z{P) = Tr{e-^^} = J D[U] e-^(^) , (11) 

where we suppress the dependence of Z on a and D[U] 
is the Haar measure over the gauge links Ufi{xo,x), for 
/i = 0, 1, 2, 3 and for each lattice point (a;o, x). For def- 
initeness we choose the action S{U) to be the standard 
Wilson plaquette action , whose bare coupling we de- 
note by (for unexplained notation see Ref. |10|). 

The most straightforward way to compute 
Z{/3,z)/Z{f3) as given in Eq. ^ is to define a set 
of (n + 1) systems with partition functions Z{/3,ri), 
(ri — i/n, i = 0,1,..., n) and corresponding actions 
S{U,ri) = r,S{U) + (1 - r,)5([/^), where is the field 
obeying shifted boundary conditions, see Eq. ([5]). We 
can then write the identity 



Continuum extrapolation of s/T^ 
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FIG. 1. Scaling behavior of s/T^, see Eq. The Stefan- 

Boltzmann value reached in the high-T limit is also displayed. 

If one defines the 'reweighting' observable as 

0(C/,r,+i) = e^(^^'-'-^'-^(^-'-'', (13) 

then the ratio Z{f3, ri)/Z{(3, r^+i) in Eq. (fT2|) can be com- 
puted as its expectation value on the ensemble of gauge 
configurations generated with the action S{U,ri+i) |lO| . 
The generating functional can thus be written as 

-(ft-)-|">.{^}- (") 

The continuum limit value of the second cumulant is ob- 
tained from a discrete lattice by a limiting procedure, 

where z = (71^0,0,0), the integer Uz being kept fixed 
when a 0. The continuum value of the entropy is thus 
approached with O(a^) corrections. 

If the set of interpolating systems is chosen so that 
the error of each contribution in the sum is comparable 
and n — {L/a)^, then the cost of the simulations for a 
given relative precision on K{f3, z, a), at fixed value of z, 
scales approximatively as a~^, while it remains roughly 
constant as a function of L. Each derivative of K at the 
origin requires an extra factor oc in statistics and 
therefore in the overall cost. Such a cost figure does not 
take into account autocorrelation times, and it depends 
heavily on the particular algorithm that we have imple- 
mented. For instance, a more refined interpolation path 
between the initial and final systems may lead to a sig- 
nificant speedup. 

We have calculated the entropy at three temperatures, 
1.5, 4.1 and 9.2Tc, see Table U for the numerical results. 
The update algorithm used is the standard combination 
of heatbath and overrelaxation sweeps [l9l - l2]| . The only 
changes over the standard algorithm reside (a) in the 
computation of the 'staples' that determine the contri- 
bution to the action of a given link variable U^j,{x) and 
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(b) in the more frequent updating of the two time-shces 
on which the observable has its support. 

The bare coupling was tuned using the data of 
in order to match lattices of different /3/a to the same 
temperature. Motivated by a study of the free case, we 
chose nz — 2. The scaling behavior of the entropy den- 
sity is displayed in Fig. [T] We observe that the cutoff ef- 
fects are quite mild. Taking into account the systematic 
uncertainty in performing the continuum extrapolation, 
our results at the two lower temperatures are compatible 
with the published results [H, [23|. The results at 9.2Tc 
are new. 



V. FINAL REMARKS 

In this Letter we have introduced the generating func- 
tion of cumulants of the total-momentum distribution 
and a new way of computing it. As an application, we 
have shown that the second cumulant is a measure of the 
enthalpy (or particle density in the non-relativistic case, 
see Eq. ([8])), and calculated the enthalpy density of glu- 
ons at three different temperatures. The ideas presented 
in this paper have further interesting applications. 

One application is the determination of the partial 
pressures of different symmetry sectors (labeled by con- 
tinuous as well as discrete quantum numbers) in the con- 
fined phase of QCD. This will provide a far more strin- 
gent test of the hadron resonance gas model than has 
been possible so far. The latter model postulates that the 
QCD pressure is due to the sum of the partial pressures 
of all zero-temperature resonances of width F < T, and is 
an important ingredient in the phenomenology of heavy- 



ion collisions (see for instance [25|). Beyond the pressure, 
any observable can be calculated in the restricted ensem- 
ble where certain quantum numbers assume prescribed 
values. 

Our method has a 'kinematic' character, and the lat- 
tice action is the only ingredient in the calculation. It is 
conceivable that the boundary conditions adopted in this 
paper can be used to formulate Symanzik improvement 
conditions |9j , or to compute the constants needed to 
define an energy-momentum tensor which satisfies the 
Ward identities in the continuum limit [2^. 

We stress that the generating function (/?, /i, z) is of 
intrinsic interest, beyond giving access to basic thermo- 
dynamic quantities. It can be used to assess how nearly 
scale-invariant the system is at a given temperature, see 
Eq. pn]) . In QCD at finite baryon density and in the 
low-temperature limit, the generating function K is an 
order parameter for the spontaneous breaking of transla- 
tion invariance. Finally, it may be of interest in analytic 
treatments, where it is more convenient to absorb the 
shift z into the action [U. The propagators are then 
modified, the parameter iz/ P playing the role of an ex- 
ternal velocity parameter. 
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